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A simple model of an atomic Bose-Einstein condensate in a box whose size varies with time 
is studied to determine the nature of adiabaticity in the nonlinear dynamics obtained within the 
Gross-Pitaevskii equation (the nonlinear Schrodinger equation). Analytical and numerical methods 
are used to determine the nature of adiabaticity in this nonlinear quantum system. Criteria for 
validity of an adiabatic approximation are formulated. 
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I. INTRODUCTION 

The Adiabatic Theorem of quantum mechanics insures that an eigenstate of a system whose Hamiltonian evolves 
sufficiently slowly in time (as determined by criteria for the applicability of the theorem) will remain in the same 
eigenstate, even though the eigenstate evolves in time Hence, a slowly evolving system which is initially in its 

ground state will remain in the ground state throughout the course of its evolution. The adiabatic theorem relies 
heavily on the superposition principle of quantum mechanics (although in classical mechanics similar theorems are valid 
for nonlinear systems [Q ) . It is of interest to determine to what extent adiabaticity carries over to nonlinear quantum 
systems, such as Bose-Einstein condensates (BECs) in the region where the mean-field description is appropriate. 
Well below the critical temperature, the mean-field description is based on the Gross-Pitaevskii equation (GPE), 
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and this approximation often yields excellent results for the system dynamics, even when the external potential V 
varies with time. In Eq. (|l|), = ATTaofi^ /m is the atom-atom interaction strength that is proportional to the s-wave 
scattering length oq, and m is the atomic mass. The parameter A'o in Eq. (|l|) is the total number of atoms, and the 
wave function 4* is subject to the normalization J |^(r, t)\'^dr = 1 (the normalization integral is a dynamical invariant 
of GPE). Adiabatic considerations regarding the GPE dynamics have been applied to cold Bosonic atoms trapped 
in optical lattices 1^-^, and the formation of optical lattice gates for quantum computing from atomic BECs 
However, the applicability of the adiabaticity concept to BECs does not follow from the above-mentioned Adiabatic 
Theorem of quantum mechanics, since the nonlinearity does not allow applicability of the superposition principle to 
the GPE. On the other hand, adiabaticity of nonlinear wave equations, and in particular, of soliton solutions to such 
equations, have been extensively studied (for a review, see Ref. jl^). Nevertheless, BEC dynamical problems based 
on the GPE have their own specific features, so that this case can be different from that studied in the framework 
of perturbed soliton solutions to the nonlinear Schrodinger equation (NLSE) in other contexts. Here, we develop 
a physically relevant one-dimensional BEC model which we study in detail by means of analytical and numerical 
methods to determine the nature of adiabaticity in nonlinear quantum systems within mean-field theory. 

There are several regimes in which adiabaticity can be experimentally and theoretically probed for nonlinear systems. 
The simplest regime is one for which the characteristic dynamical time scale (i.e., the time during which parameters 
of the Hamiltonian undergo an essential change), T, satisfies conditions 



TAD < r < TNL- 



(2) 



Here, tad is the quantum-mechanical linear adiabatic time scale determined in terms of the inverse of the difference 
of the energy eigenvalues at different values of time, tad = max{?i/ [ei{t) — eo(t)]}, where the maximum is taken 
with respect to a given time interval while the nonlinear time scale is tjvl = nia,x{h / ^{t)} , with /i(i) being the 
instantaneous chemical potential jllj. In this case, the applicability of the linear Adiabatic Theorem is secured by 
the first inequality in (0) , and nonlinearity cannot play a significant role in the dynamics due to the second inequality. 
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Therefore, the dynamics must be adiabatic. In particular, this regime apphes to the NIST optical-lattice experiments 
wherein light pulses in the jis range are applied to a sodium BEC |l2) . 

A more interesting and more problematic regime is when the dynamical time scale is the largest one in the problem, 
i.e., tadjTnl T. This case applies, e.g., to the BEC experiments reported in Refs. ^,|^. Here, the nonlinearity 
plays an essential role in the dynamics, and a relevant question is whether the dynamics can be adiabatic. Generally, 
a further time scale appears in the multidimensional GPE, viz., the diffraction time, tjjf — 2mL^/h, where Lq is 
the length of the system (see below). In what follows, we show that the GPE does allow for adiabaticity when 
TAD, tdp^tnl ^ T, and we give explicit criteria for the validity of adiabaticity in the GPE. We do so by means of an 
analytical estimate of corrections to the adiabatic approximation, and we present numerical results for the dynamics 
of this system. 



II. THE DYNAMICAL MODEL 



We consider a model based on the ID GPE, in which the external potential U{x) is an infinitely deep well, 

[/(.)= '^<^'^^'^m (3) 

^ ^ [ +00, X <0 01 X > L{t) . ^ ' 

The size of the well L[t) slowly varies with time, and we are interested in determining the behavior of the system in 
this case, to determine the applicability of adiabaticity. The ID GPE takes the form 

fi^ 2 

ihipt^ -i^'^xx + g\M tp , (4) 

with the boundary conditions 

V(0,t) =V(i(t),i) =0 , (5) 

and normalization of the wave function, 

\ij{x)\^dx = l. (6) 

The nonlinearity parameter g appearing in this ID GPE is related to the nonlinearity parameter NqUq in its 3D 
counterpart, Eq. (l]), and is determined so that iVoC/ol^'mP ~ .'7|V'mP, where ^I'm and 4'm are the maximum values of 
the 3D and ID wave functions, respectively. This condition insures that the time scales for the nonlinear interaction 
in the 3D and ID cases are equal (see below and Rcf. 

This is the generalization of the particle in a box problem to the case where (a) the size of the box is varying with 
time, and (b) there are many bosonic particles in the box that arc interacting via a mean field. 

A typical situation in which the dynamics may be adiabatic is when the function L(t) takes on constant values as 
t — > ±00, and slowly varies in between on a long time scale T. We aim to find the final state tpix, t — +00) into which 
an initial state 'ip{x,t = —00) will be transformed if T is sufficiently large, and to check whether the wave function 
ij}{x,t) remains adiabatic during the course of the evolution, provided that the function L(t) varies slowly enough (in 
practice, of course, the evolution time interval is large but finite). To determine what "sufficiently slow" means, we 
define the nonlinear time scale obtained directly from the GPE as 

tNL = {9\i'^?in)-^ « ^i/n , (7) 

where /i and |'0m| are the chemical potential and maximum of the wave function in the initial configuration. The 
evolution is is slow as compared to nonlinear time scale if T ^ t^L- In many BEC systems, the nonlinear time scale 
is large compared with the diffraction time scale tuf = 2mLo/7i, also obtained directly from the GPE so we 
should have T ^ tj^L ^ tdf- 

It is convenient to transform the variables x, t and ip to new (dimensionless) variables r, ^ and u: 

e = ^/m , (8) 

r=^r^ (9) 
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(10) 



Note that the problem is mapped onto a fixed spatial interval ^ G [0, 1] of the dimcnsionless spatial variable f , and 
the boundary condition is therefore not time dependent when the problem is reformulated in terms of these variables. 
The ID GPE (^ takes the following form in terms of the new variables: 

iu-T + U£^£^ — \u\^ u = i{Lr/L){S^u)^ , (11) 

where L-r = dL/dr. Equation (11) is supplemented by boundary conditions following from Eq. (^, 

u(e = 0,r)=u(^=l,r) = 0. (12) 

The norm defined in terms of the transformed wave function u, 

N[u{t)]= ['\u{^,rfd^, (13) 

is not conserved in time, unlike the original norm, J^^*^ I'fpi^i ^)l^ d,x — 1. Indeed, as follows from the substitution of 
Eq. ( p^ ) for u into Eq. (^3|), the u-norm is an explicit function of time: 



N[u{r)]='-^Lir)^no^, (14) 



where Lg — L{t = —oo) (or alternatively Lq = L(t — tg) if the initial moment in time is Iq). The dimcnsionless 
nonlinear-strength parameter, 

no = 2gmLo/h^ , (15) 

introduced in Eq. (^) will play an important role below. 

When the system size L(t) is a slowly varying function of time, the right-hand side (RHS) of Eq. dll] ) is small, 
being proportional to the logarithmic derivative of the slowly varying function. Therefore Eq. ( [ll| ) may be naturally 
considered as a perturbed self-defocusing NLSE, and the adiabatic methods for nonlinear wave equations reviewed in 
Ref. jl^ might be applied. However, the perturbation term on the RHS of Eq. ( [ll| ) need not allow straightforward 
application of the perturbation theory to the present problem since this term does not vanish at ^ = and ^ = 1 
when a general solution found in the zeroth-order approximation (the expression ( p^ ) below) is inserted into it. One 
can easily check that, as a consequence of this, a perturbative expansion generated by the term on RHS of Eq. (O) 



is incompatible with the boundary conditions (12) 



To resolve the problem, we transform the wave function once again, defining 

2 L-^ ^2 



u{^,T)^vitT) e^p{-^e]- (16) 



The transformation (|l^) generates a more convenient form of the perturbed NLS equation 

,2 -Lt LLtt ~ 2i^ 2 



.Lr Z/Lt-t — 21/^ 2 

v\ V = i—v H C " I (1 ') 



2L 4L2 



which is subject to the same boundary conditions as in Eq. (|12|) , = 0, t) = = 1, t) = 0. An obvious advantage 
of having the perturbed NLS equation in the form (|l^) is that now the perturbation vanishes at ^ = and ^ = 1, 
once a solution found in the zeroth-order approximation vanishes at these points. 

Note that the first term on the right-hand side of Eq. (^7|) is non-conservative. Accordingly, it is straightforward to 
see that this term leads to the exact relation ( |l^ ) for the norm evolution. Another important fact is that the second 
term on the RHS of Eq. ([T^), unlike the first term, is second-order small with regard to derivatives of the slowly 
varying functions. In the perturbation-theory section that follows below, we will not consider effects produced by the 
second-order term, focusing solely on the most important first-order effects. 
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III. ADIABATIC PERTURBATION THEORY 



A. Zeroth-Order Approximation 

In zeroth-order approximation of the perturbation theory (neglecting the RHS of Eq. (|l7|)), an exact stationary 
solution satisfying the zero boundary conditions at ^ = and ^ = 1 is given by |]l3|-|l5|: 

v{^, t) = 2^l'^kK{k) sn{2K{q)^, k) exp(-i^r) = V{^] k) exp {icj) (r)) . (18) 

Here sn(-, •) is the doubly periodic Jacobi elliptic sine function, k is the corresponding elliptic modulus, K{k) is the 
complete elliptic integral of the first kind, and the chemical potential is related to k as follows: 

^ = A{l + k^)K^{k). (19) 

The modulus fc, which takes values < fc < 1, determines the strength of the nonlinearity: it is weak if fc — > 0, 
and strong if fc 1. In fact, k is related directly to the dimensionless nonlinearity-strcngth parameter no defined in 
Eq. ( p^ ) as follows: 

8K{k) [K{k) - E{k)] = no , (20) 

where E{k) is the complete elliptic integral of the second kind. Thus, k completely determines the normalization of 
the initial wave function, and visa versa. 

To illustrate the zeroth-order solution, plots of 8K(k) [K{k) — E{k)] versus k (see Eq. (pO|)), and 
2'^^^kK{k) sn(2_ft'(g)^, k) versus ^ for three different values of k in the regime of strong nonlinearity, k > 0.5 (see 
Eq. ([l^)), are displayed in Figs. |l|(a) and 0(b). We remark that an exact solution that can be expressed in terms of 
the Jacobi elliptic functions is frequently called a cnoidal wave, which stems from the notation cn for the Jacobi's 
elliptic cosine, related to the elliptic sine. 



B. The Nonlinear Adiabatic Approximation 

The first-order perturbation term on the RHS of Eq. ( p^ can be treated in terms of nonlinear adiabatic perturbation 
theory ||l^]. We stress that, unlike the perturbation term in the intermediate equation (|Tl|), which is "abnormal" in 
the sense that it is not compatible with the necessary boundary conditions, as it was explained above, the "normal" 
perturbation in Eq. (|l^) satisfies the boundary conditions. The applicability of simple perturbative techniques for this 
class of models can be proved using a rigorous expansion based on the inverse scattering transform for the unperturbed 
NLS equation (i.e., one can prove that the "simple techniques" yield, in the lowest-order nontrivial approximation. 



exactly the same results as the rigorous methods, see Ref. |10| and references therein) 



The first standard step of the perturbative analysis is to apply the lowest-order adiabatic approximation. This 



approximation takes the unperturbed solution (18), which contains the parameter fc, and makes it the first-order 
approximate solution to the perturbed equation, assuming that the modulus k is slowly varying in time, rather than 
remaining constant. 

The slow dependence of the parameter(s) is introduced so as to cancel the secular divergence(s) in the perturba- 
tion theory. An important case is when the unperturbed solution contains a single nontrivial parameter (fc, in the 
present case), and the perturbed equation gives rise to an exact relation replacing a conservation law existing in the 
unperturbed version of the equation (this exact relation is usually called a balance equation for the (former) conserved 
quantity). This is the case in Eq. (|l4|). Then, the time dependence of the parameter, i.e., k{T), can be found in a 
very simple way by substitution of the zeroth-order approximation for the solution into the balance equation |p^ . In 
the present case, this condition amounts to evaluation of the actual value of the norm (|l3|), inserting the solution (|l^ ) 
into it, and then substituting the result into the exact relation (llj). The final result is 



8K{k) [Kik) - E{k)] = no L{t)/Lo . (21) 



Eq. (21) is a transcendental equation to determine fc(r) for a given function L(t) and no (recall that no is a constant). 

An essential ingredient of the adiabatic approximation is a consistent definition of the phase 0(r) for the first-order 
solution with variable A:(t). Indeed, substituting fc(r) back into the general expressions ( |l8| ) and ( p^ ) for the wave 
function, it is easy to see that the consistently defined phase is 
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0(t) = - f\{T')dT' = -4 r [1 + e{T')] (fc(T')) dr' , 



(22) 



To being the initial time (tq = — oo in the usual formulation of the adiabatic approximation). 

Thus, the full expression for the lowest-order perturbative solution obtained in the adiabatic approximation is 

v[^,T)^V{tk{T)) exp(*(/.(T)) , (23) 

where F(^; k) and (^{t) are given by Eqs. (|l^) and (p2|), respectively. Note that expression (^3|) automatically satisfies 
the zero boundary conditions at the points ^ = and ^ = 1. 

Knowing a particular form of the slow temporal dependence fc(r) obtained from Eq. (^l|), one can find the temporal 
dependence of the solution's amplitude, 

^(r) EE max|i;(C,T)| = 22/3A;(r)i^ (fc(T)) . (24) 

The temporal dependence of state's width (which, for instance, can be defined as the full width at half-maximum of 
|i^,r)p) can similarly be obtained in the adiabatic approximation from the above expressions. Using Eqs. ( |l9|) and 
(pi|), it is also possible to predict the evolution of the instantaneous value of the chemical potential /i(r). 



C. Corrections to the Lowest-Order Adiabatic Approximation 



Once the slow time dependence of fc(r) has been determined as described above, one can look for perturbation- 
induced corrections to the state's shape, which was not taken into account in the first-order adiabatic approximation. 
A solution to Eq. ( |l7| ) including the corrections can be sought in the form of an expansion compatible with the zero 
boundary conditions, namely. 



bmir) sin(7rm^) 



exp {i(j>{T)) 



(25) 



where the functions V and (fi are those which were obtained in the previous subsection. 

The simplest way to derive evolution equations for the amplitudes bm(t) is to directly substitute the expansion 
( p5| ) into Eq. (|l^), multiply the resulting equation by sin(7rm^), and integrate from ^ = to ^ = 1, carrying out this 
procedure for each integer m. The correction terms are neglected when substituting the expression ( p5| ) into the first 
perturbation term on the RHS of Eq. (^, as they would give rise to higher-order perturbations. Implementing this 
procedure, we use the classical Fourier expansion for the function sn, 



sn(2X(/c)C,fc) 



27r 



kK{k) 1 - Q^P 
Q{k) = exp 



oo QV-l/2 



t:K (VI - k"^) 



K{k) 



(26) 
(27) 



A complicated system of inhomogeneous linear evolution equations for bm{T) ensues. If m is odd, i.e., m = 2p — 1, 
we obtain 



(28) 



and, if m is even, i.e., m = 2p, 



(29) 



db 

«^ + ^2p62p - E ^'^2p,„ (26„ - &;) = . 

n=l 

Here Q is the Jacobi parameter defined in Eq. (p7|), and the coefficients appearing on the left-hand sides of Eqs. (|2. 
and (|29|) are 



5 



Mmn = Ak^K^{k) I an^{2K{k)£_, k) [cos (7r(m - n)) - cos (7r(m + n))] 



(30) 
(31) 



In fact, all the coefficients Mmn with m and n having opposite parities are zero, hence we may set b2p = 0, and we are 
left with the system of equations (p8|). Recall that one should substitute the time-dependent modulus A:(t) as found 
from Eqs. (21)) into the above expressions. 

This cumbersome system can be simplified if fc(r) does not take on values too close to unity. Then, sn remains close 
to the usual sine function (for instance, at fc^ = 1/2, the Jacobi parameter, which determines the anharmonicity of 
the expansion (26), is Q = exp(— tt) « 0.043, which may be regarded as a sufficiently small expansion factor). Thus, 
to obtain a simple approximation for the coefficients Mmn defined in Eq. (^l|), one may simply set sn(2i^(fc)^, /c) 
sin(7r^). Within this approximation the only nonzero components of in the matrix (Mmn) are 



Mil = 3k^K^{k), M,„„ = 2k^K^{k) {m > 1), M, 



= M, 



7n — 2. in 



= ^k^K^{k) 



(32) 



Furthermore, the RHS of Eq. (J^^ also greatly simplifies in the same approximation. It is different from zero solely 
for m = 1, being equal to iLr/{2L). Thus, the approximation which replaces the sn function by the usual sine leads 
to the following equations, instead of Eqs. (28) and (p9|): 



dbi 



2 (1 - 2k') K\k) - 



bi + 3eK\k)bl + k'K' (263 - &3) = ^ , 



(33) 



'2p-l 



dr 



2{l-k')K'{k)- 



W2p-i)F 



2k'K'{k)b;^^ 



k^K' {2b 



'2p-3 



262P+1 — b\ 



— h* 
-3 02p4 



= 0, 



(34) 



where p > 1. Recall that all the amplitudes bm with even values of m are zero. 



Despite the fact that the approximate system consisting of Eqs. (33) and ( p4D is considerably simpler than the exact 
Eqs. (|2^), it can only be solved numerically by truncating the system of the linear equations at some finite integer. 
Nevertheless, some qualitative generic features of the solution can be determined. The general structure of the system 
is of the form 



dB 

i— + A{T)B = iC{T) , 
dr 



(35) 



where B is a column vector of the variables 62^-1, A is a matrix of coefficients multiplying the variables 6, and C is 
the vector column of free terms on the left-hand side, with the single nonzero entry ci = Lr/{2L). Both C and A 
slowly depend upon time - the former directly, the latter via fc(r). 
Solutions to the system (|3q) consist of terms of the type 



exp 



u;{r")dr" f{r')dT' , 



(36) 



where w(t) is an eigenvalue of the matrix A(t), and /(r) are slowly varying functions similar to the above-mentioned 
ci. Note that the time scales 2tt/uj, determined by different eigenfrequencies co, arc, in fact, a mixture of the adiabatic 
and nonlinear time scales, tad and tjvl, defined in the Introduction. 

The following conclusion can be made concerning the size of the nonadiabatic effects (shape corrections) considered 
above. If the function L{t) slowly depends on r with a characteristic time scale T (as defined in the introduction), 
and if a characteristic value of lu is luq (within the limits of its slow evolution on the time scale ^ T), the criterion 
for the applicability of the adiabatic approximation is 



ujoT/{2Tr) > 1 



(37) 



We stress that, as the characteristic times 2tt/ujq taken for the different eigenfrequencies constitute a set including 
the adiabatic and nonlinear time scales tad and tnl (see above), the inequality ( p7| ) is exactly the condition for the 
applicability of the adiabatic approximation conjectured in the Introduction. 
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The evolution equations ( JSq ) for the shape-correction amphtudes are to be solved for an initial state without shape 
corrections, i.e., bmir = tq) = V m. If one takes the initial moment as tq —> —oo (as mentioned above, this is the 
standard assumption in the treatment of adiabatic processes [QQ), one can determine eventual values of the shape- 
correction amplitudes as &m(T) at r — > +00. Classical estimates for integrals involving products of rapidly and slowly 
varying functions ^ show that the values of bm{T — *■ -l-oo) are exponentially small when condition ( p7| ) is satisfied: 

\b,n{T = +(X))\ ^ exp {—const ■ ujoT) . (38) 

A particular value of the constant in this expression depends on the choice of the unperturbed state and on the form 
of the function L{t). Hence, in analogy with the well-known theorems estimating nonadiabatic corrections to the 
adiabatic approximation in (nonlinear) classical mechanics the 6m (t — > -|-oo) values are exponentially small. 



IV. NUMERICAL RESULTS 



We first present results for the amplitudes 6to(t) in Eq. < \25\) , obtained by numerically solving Eqs. { p3\ ) and (|3; 
We take P = 0.3, Lq = 1, and 



Ht) = Lo 



1 + exp 



(r-T/2) 



2a2 

with a = T/IO. Figure |2| shows the computed excited-state probabihty, 



P,.(r)^^|6m(r)p, (39) 



m=l 



versus time r, with the number of modes kept in the truncated calculation being Nc = 1, 3, 5, 7, and 9, for T = 100. 
Except for A^c = 1, all the curves lie on top of each other, hence the results do converge very quickly as a function of the 
number of the modes. We see from Fig. ^that the probability of finding excited states for all times is below 3.2 x 10^^, 
and for t = T the probability is exceedingly small, i.e., the process is almost completely adiabatic. A minimum of 
Pcx{t = T/2) is expected from the general form of the perturbation equations since the derivative of L{t) vanishes at 
T = T/2. For T ^ 10, the excited-state probability ( |39| ) begins to be large (> 0.2), and the adiabatic-theory results 
are no more reliable. For example. Fig. ^ shows the results for T = 10. Again the convergence as a function of the 
number of the modes is very fast, but the excited state probability is not small. For times r > T, Pcx{t) oscillates 
with time. 

For stronger nonlinearity, 0.5 < /c^ < 1, perturbative methods cannot be used (in particular, the approximation 
based on the replacement of the elliptic sine by the ordinary sine, as described in the above section, does not apply), 
so we directly solved Eq. ( p7| ) using a split-step fast Fourier transform method, in order to check if adiabaticity still 
takes place in this regime. Figure ^ shows the results for the calculated wave function ip{^,T = T) = \ip{^)\e^^'^^^ 
versus ^ in the box at the completion of the dynamical process for small (non-adiabatic) time-scale, T = 0.01, and 
for k — 0.963. Also shown for comparison is the initial eigenstate magnitude \^|J{^,T = 0)|. The magnitude of the 
wave function in the final state is not too different from that in the initial eigenstate, and the spatial variation of the 
phase is fairly flat. Figure js] pertains to the same case, but with a larger time scale, T = 1. Now, the magnitudes of 
the wave function in the finite and initial state are barely distinguishable on the scale of the figure, and the spatial 
profile of the phase is almost flat. Thus, the process is largely adiabatic in the latter case. 



V. CONCLUSION 



We have presented a consistent derivation of the nonlinear dynamics for a simple model describing a BEC confined 
in a box with a temporally- varying size L{t). The speed of the variation of L(t) determines whether the dynamics 
is adiabatic. A "trivial" regime of adiabaticity is that for which tad T t^vl! in this work, we have shown 
that adiabaticity can also be maintained when tad,tdf,tnl <C T, where the various time scales have been defined 
in the Introduction. If other time scales appear, the condition tad,tdf,tnl ^ T may not be sufficient to insure 
adiabaticity. For example, if a barrier is present in the middle of the box - e.g., a repulsive delta- function at x = Lq/2 
- then another time scale, corresponding to the time of tunneling under the barrier, r^, is present in the problem. If 
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this time scale is long compared with tjvl , adiabaticity will not be maintained, and a non- vanishing spatially- varying 
phase will develop across the condensate wave function [ p^ . Hence, the issue of adiabaticity in nonlinear problems 
must be investigated carefully; the perturbative techniques reviewed in Ref. may be applicable, but additional 
considerations may play a role. 

The particle in a box is a paradigm problem in one-body quantum mechanics. We have extended it to the many- 
body regime at least within a mean-field approach, and studied the adiabaticity for such a system when the size of 
the box varies with time. Specifically, we formulated the criteria for the validity of the adiabatic approximation for a 
BEC in a box whose size varies with time, developed the analytical and numerical tools for investigating adiabaticity 
in the dynamics within quantum mean-field theory, and presented results of calculations for this system. 
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FIG. 1. (a) The expression 8K{k) [K{k) — E{k)] versus the elhptic modulus k. (b) The zeroth-order analjrtic solution 
2^^^kK(k) sn(2K{q)^, k) for three diflterent values of k. The normalization of these soliton solutions can be read off the curve 
in (a). 
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FIG. 2. The probability X/mLi foi' the nonadiabatic correction to the state versus r for Nc = 1,3,5,7,9 with 

T = 100. For A^c > 3 the curves lie on top of each other. 
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FIG. 5. The magnitude and phase of the wave function ip{(,,T — T) — versus the coordinate £, in the box at the 

completion of the the dynamical process, with k — 0.963 and T = 1. 
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